Black Hole Excision for Dynamic Black Holes 
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We extend previous work on 3D black hole excision to the case of distorted black holes, with a 
variety of dynamic gauge conditions that respond naturally to the spacetime dynamics. We show 
that in evolutions of highly distorted, rotating black holes, the combination of excision and the gauge 
conditions we use is able to drive the coordinates to a frame in which the system looks almost static 
at late times. Further, we show for the first time that one can extract accurate waveforms from these 
simulations, with the full machinery of excision and dynamic gauge conditions. The evolutions can 
be carried out for long times, far exceeding the longevity and accuracy of better resolved 2D codes. 
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Introduction. The long term simulation of black hole 
(BH) systems is one of the most challenging and impor- 
tant problems in numerical relativity. For BHs, the dif- 
ficulties of accuracy and stability in solving Einstein's 
equations numerically are exacerbated by the special 
problems posed by spacetimes containing singularities. 
At a singularity, geometric quantities become infinite and 
cannot be handled easily by a computer. 

Traditionally, in the 3+1 approach the freedom in 
choosing the slicing has been used to slow down the ap- 
proach of the time slices towards the singularity ( "singu- 
larity avoidance"), while allowing them to proceed out- 
side the BH . Singularity avoiding slicings are able to 
provide accurate evolutions, allowing one to study BH 
collisions and extract waveforms but only for limited 
evolution times. Combining short full numerical evolu- 
tions with perturbation methods, one can even study the 
plunge from the last stable orbit of two BHs But a 
breakthrough is required to push numerical simulations 
far enough to study orbiting BHs, requiring accurate evo- 
lutions exceeding time scales of t w 100M. In 3D, tradi- 
tional approaches have not been able to reach such time 
scales, even in the case of Schwarzschild BHs. Character- 
istic evolution codes, on the other hand, are well-adapted 
to the long-term evolution of single black holes M, but 
here we concentrate on non-characteristic methods for 
their ease of generalization to multiple black holes. 

A more promising approach involves cutting away 
the singularity from the calculation ("singularity exci- 
sion"), assuming it is hidden inside an apparent hori- 
zon (AH) (|||. This work has been progressing, from 
early spherical proof of principle || to recent 3D de- 
velopments |0,^-^O[. However, beyond a few spherical 
test cases [ pT[p2| , excision has yet to be used in conjunc- 
tion with live gauge conditions designed to respond to 
both the dynamics of the BH and the coordinate motion 
through the spacetime. 

In this paper we extend recent excision work [|| to the 
case of distorted, dynamic BHs in 3D, using a new class 
of gauge conditions. These gauge conditions not only re- 



spond naturally to the true spacetime dynamics, but also 
drive the coordinates to a frame where the system looks 
almost static at late times. We show that not only are the 
evolutions accurate as indicated by the mass associated 
with the apparent horizon, but also that very accurate 
waveforms can be extracted, even when the waves carry 
only a tiny fraction of the energy of the spacetime. We 
also show that the 3D evolutions of dynamic BHs we are 
now able to perform, are superior, in terms of accuracy, 
stability, and longevity, to previous 3+1 BH simulations, 
whether carried out in full 3D or even when restricted to 
2D. These results indicate that BH excision can be made 
to work under rather general circumstances, and can sig- 
nificantly improve the length of the evolutions, and the 
accuracy of the waveforms extracted, which will be cru- 
cial for gravitational wave astronomy. 

Initial Data. For thi s pap er we consider single dis- 



torted BH spacetimes 13,14] that have been used to 
model the late stages of BH coalescence (j^pl. Fol- 
lowing |TJ,|lJ], the initial 3- metric j a b is chosen to be 
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where the "Brill wave" function q is a general function of 
the spatial coordinates, subject to certain regularity and 
fall off restrictions, that can be tailored to produce very 
distorted 3D BHs interacting with nonlinear waves. The 
radial coordinate r\ is logarithmic in the cartesian radius 
r. There are two classes of data sets used here corre- 
sponding to even- and odd-parity distortions. The even- 
parity data have vanishing extrinsic curvature, while the 
cases containing an odd-parity component have nontriv- 
ial extrinsic curvature Kij. As shown in |l7| , |l8|] , these 
distorted BH data sets can include rotation as well, corre- 
sponding to spinning, distorted BHs that mimic the early 
merger of two orbiting BHs. Hence they make an ideal 
test case for the development of our techniques. We leave 
the details of the construction of these BH initial data 
An important point that we wish 



sets to Refs. 17 



to emphasize is that such data are not of the Kerr-Schild 
form with ingoing coordinates at the horizon. That form 
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of initial data sets has been recently advocated since it 
is not conformally flat Jis[| and is well adapted to inward 
propagation of quantities at the horizon. 

Evolution and Excision Procedures. Our simulations 
have been performed using what we refer to as the 
"BSSN" version of the 3+1 evolution equations [p)|-p||, 
which we have found to have superior stability properties 
when compared to standard formulations. As detailed 
in p2|-p4||, we actively force the trace of the conformal- 
traceless extrinsic curvature Aij to remain zero, and 
we use the independently evolved "conformal connec- 
tion functions" T l only in terms where derivatives of 
these functions appear. All the simulations described be- 
low have been performed using a 3-step iterative Crank- 
Nicholson scheme and a radiative (Sommerfeld) outer 
boundary condition. We refer the reader to Ref. |24| for 
the details of the numerical implementation. 

We use the simple excision approach described in 
Our algorithm is based on the following ideas: (a) Excise 
a cube contained inside the AH that is well adapted to 
cartesian coordinates; (b) Use a simple boundary con- 
dition at the sides of the excised cube: copying of time 
derivatives from their values one grid point out along 
the normal directions; (c) Use centered (non-causal) dif- 
ferences in all terms except for advection terms on the 
shift (terms of the form P l di ). For these terms we use 
second order upwind along the shift direction. These 
simplifications in excision reduce the complexity in the 
algorithm, avoid delicate interpolation issues near the 
excision boundary, and have allowed us to make rapid 
progress. Currently, the method is implemented for non- 
moving excision regions, although they are allowed to 
grow. One can hope that colliding black holes can be 
treated even with this restriction through the use of co- 
moving coordinates. A more detailed description of our 
excision algorithm can be found in Ref. 

Gauge Conditions. For the lapse we use a hyperbolic 
slicing condition motivated by the Bona-Masso family of 
slicing conditions po| , 

d t a = ~a 2 f(a) (K - K ), (2) 

where K is the trace of the extrinsic curvature, K$ is its 
value in the initial data, and / is a (positive) function 
of a which we specify below. With this condition, the 
lapse will evolve as long as a 2 f(a) and K — Kq are non- 
vanishing. 

For the shift f3 z we have considered families of elliptic, 
parabolic, and hyperbolic conditions that relate the shift 
with the evolution of the conformal connection functions 
r*. We obtain parabolic and hyperbolic shift conditions 
by making either dt(3 l or df(3 l proportional to the ellip- 
tic operator for (3 l contained in the "Gamma freezing" 
condition d t f k — (see Ref. f24|l ), itself closely related 
to the well known minimal distortion family Q. Ellip- 
tic conditions have the disadvantage of requiring bound- 
ary data at the excision region where it is difficult to 
know what to impose, while parabolic conditions force 



a strong restriction on the stability of the differencing 
scheme: At oc (Aa;) 2 (this is true for explicit schemes, 
implicit schemes have no such restriction) . We have then 
concentrated on hyperbolic conditions of the form 

d?p = cd t r-td t p\ (3) 

where C and £ are positive functions. We call such evo- 
lution conditions for the shift hyperbolic "driver" condi- 



In the spirit of the puncture method for evolutions 




we use a BSSN scheme with the usual time-dependent 
conformal factor e 4 * and an additional time-independent 
conformal factor ^ 4 that comes from the initial data. In 
our examples we use £ = fc/'I' 4 , where k is a positive 
constant. The division by ^ 4 helps to slow down the 
evolution of the shift in the vicinity of the black hole. We 
have found it important to add a dissipation term with 
a constant coefficient £ in order to reduce some initial 
oscillations in the shift. Notice that in contrast with k, 
the coefficient £ is not dimensionless (it has dimensions 
of inverse length), so in practice we rescale it using the 
total mass of the system. Experience has shown that by 
tuning the value of £ we can almost freeze the evolution 
of the system at late times. 

The parameters used for all simulations described be- 
low are: a is given by Eq. (0), with / = 2/ a; j3 l is given 
by Eq. (§) with C = 0.75/** £ = 3/M (M is the initial 
ADM mass of the system) . As initial conditions we take 
a = 1, [3 l = 0, d t a = d t fi l = 0, except in one case men- 
tioned below where we perform a single maximal slicing 
solve to obtain a more appropriate initial lapse. Given 
these initial conditions, we let the gauge conditions take 
care of the rest. We use the same gauge parameters for 
all results in this paper, whether they are applied to 
Schwarzschild, distorted, or rotating BHs, showing the 
strength and generic nature of these conditions. 

Results. The first example we show is Schwarzschild, 
written in the standard isotropic coordinates used in 
many BH evolutions. Note that with this initial data 
and our starting gauge conditions, the BH should evolve 
rapidly. If a and (J 1 were held fixed at their initial values, 
the slices would hit the singularity at t = ttM. Instead, 
a and ft 1 work together with excision to rapidly drive the 
coordinates to a frame where the system looks essentially 
static, corresponding to the true physical situation. 

In Fig. [j] we show the radial metric function g rr /'^ i 
vs. time. The grid covers an octant with 128 3 points 
(Ax = 0.2, M — 2). Appropriate symmetry conditions 
are applied on the faces of the octant for the different 
dynamical variables. We have checked that removing 
the octant symmetry (while using a lower resolution) 
does not change the results for the evolution times re- 
ported here (in particular no instabilities were encoun- 
tered, cmp. H). Notice that the metric begins to grow, 
as it does without a shift, but as the shift builds up the 
growth slows down significantly. At this stage, the sys- 
tem is effectively static, even though we started in the 
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FIG. 1. We show the evolution of the radial metric function 
g rr /'fy i for a Schwarzschild BH along the x— axis, constructed 
from the cartesian components. The upper panel shows the 
grid stretching in the metric for singularity avoiding slicing 
with vanishing shift and no excision, while the lower panel 
shows the metric for the new gauge conditions with an exci- 
sion box inside a sphere of radius 1M. Note the difference 
in the vertical scales. Without shift and excision the metric 
grows out of control, while with shift and excision a peak be- 
gins to form initially but later freezes in as the shift drives 
the metric to a static configuration (note the time labels). 

highly dynamic isotropic coordinates. We also show the 
time development of a and (5 r in Fig. [^, which evolve 
rapidly at first but then effectively freeze, bringing the 
whole system to an almost static state by t = 10 M. The 
evolution of the metric and gauge variables then proceeds 
only very slowly with time until the simulation is stopped 
well after t = 200M. The decision to stop the code at 
this time is simply a CPU time consideration, but we 
notice that in this and all following examples the code is 
stopped once all the interesting dynamics have finished. 

Figure || shows the AH mass Mah — V 'Area^H /167r, 
determined with a 3D AH finder For comparison 

we also show the value of Mah for the 3D run without 
shift, and for a highly resolved 2D simulation with no 
shift and no excision JlJJ]. The 2D code uses maximal 
slicing, so the coordinate time t refers to different slices, 
but the slicings turn out to be even more similar than is 
to be expected from Eq. (Q). While the 3D simulation 
with shift and excision continues well beyond t — 200A/, 
the 2D result becomes inaccurate and the code crashes 
due to axis instabilities by t = 150M, and the 3D run 
without shift crashes already by t — 50M. Notice that in 
the 2D case, after around t — 35M, Mah grows rapidly 
due to numerical errors associated with grid stretching. 
With excision and our new gauge conditions, the 3D run 
has less than a few percent error by t = 200M, while 
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FIG. 2. We show the lapse and shift for the excision evolu- 
tion of a Schwarzschild BH. After around 10M, the lapse and 
shift freeze in as the metric is driven to a static configuration. 
The size of the excision box was allowed to grow with the 
change in the coordinate location of the AH. 



the 2D case has more than 100% error before it crashes 
at t f» 150Af. For the excision run, notice also that 
while there is some initial evolution in the metric and 
the coordinate size of the AH (see Figs. | and |) the AH 
mass changes very little. 

Next, we turn to a truly dynamic, even-parity distorted 
BH. This system contains a strong gravitational wave 
that distorts the BH, causing it to evolve, first nonlin- 
early, and then oscillating at its quasi-normal frequency, 
finally settling down to a static Schwarzschild BH. This 
provides a test case for our techniques with dynamic, 
evolving BH spacetimes, and allows us to test our ability 
to extract gravitational waves with excision for the first 
time. In this case, in the language of fUjl , we choose the 
Brill wave parameters to be Qo = 0.5, ?7o = 0, a = 1, 
corresponding to a highly distorted BH with M = 1.83. 
Just as before, we use a grid that covers one octant, with 
128 3 points and Aa; = 0.2 

In Fig. |] we show the AH mass Mah as a function of 
time for the distorted BH simulations carried out in both 
2D and 3D. Mah grows initially as a nonlinear burst of 
gravitational waves is absorbed by the BH, but then lev- 
els off as the BH goes into a ring-down phase towards 
Schwarzschild. Fig. || shows the proper polar circumfer- 
ence of the AH divided by its proper equatorial circum- 
ference. This ratio allows an estimate of the size of the 
local dynamics during the run. Notice how the horizon 
starts far from spherical (with a ratio close to 2), it later 
oscillates from prolate to oblate and back again, and fi- 
nally settles on a sphere (with a ratio of 1). 

In the 3D case, the gauge conditions and excision 
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FIG. 3. The solid line shows the development of the AH 
mass Mah, determined through a 3D AH finder, for the sim- 
ulation of a Schwarzschild BH shown above, while the dashed 
lines show the AH mass obtained using 2D and 3D codes with 
no shift and no excision. The 2D code crashes at t ~ 150Af, 
the 3D run without shift crashes at t ~ 50M, while the 3D run 
with shift and excision reaches an effectively static state and 
the error remains less than a few percent even after £ = 200A/. 
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FIG. 4. We show the AH masses Mah for a BH with 
even-parity distortion for the 2D (no excision, no shift) and 
3D (excision, shift) cases. The 3D result continues well 
past 150M, while the 2D result becomes very inaccurate and 
crashes by t = 100M. 



quickly drive the metric to an almost static configuration, 
as the system itself settles towards a static Schwarzschild 
BH. The evolution is terminated at around t = 160M. 
To our knowledge, distorted BHs of this type have never 
been evolved for so long, nor with such accuracy, in ei- 
ther 2D or 3D. By comparison, in the more highly re- 
solved 2D case with zero shift and no excision, the fa- 
miliar grid stretching effects allowed by the gauge choice 
lead to highly inaccurate evolutions after some time with 
the error in Mah again approaching 100% when the code 
finally crashes at t « 100M. 

In Fig. ^, we show the results of extracting waves from 
the evolution of this highly distorted, excised BH. Using 
the standard gauge-invariant waveform extraction tech- 
nique, the Zcrilli function is shown for both the 2D and 
3D simulations discussed above. There is a slight but 
physically irrelevant phase difference in the two results 
due to differences in the slicing; otherwise the results are 
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FIG. 5. We show the ratio of the polar and equatorial cir- 
cumferences as a measure of the dynamics of the AH for the 
BH with even-parity distortion . 



remarkably similar (the waves are extracted at the same 
Schwarzschild radius in both cases). This shows conclu- 
sively that the excision and live gauge conditions do not 
adversely affect the waveforms, even if they carry a small 
amount of energy (around 10~ 3 Madm in this case). 

We now turn to a rather different type of distorted 
BH, including rotation and general even and odd-parity 
distortions. In the language of Ref. jl7), the parameters 
for this simulation are Qo = 0.5, Tjo = 0, a = 1, J = 35, 
corresponding to a rotating distorted BH with M = 7.54 
and an effective rotation parameter a/M = 0.62. Pre- 
viously, such data sets could be evolved only to about 
40Af |16|. For the purposes of this paper we have chosen 
an axisymmctric case so that we can compare the results 
to those obtained with a 2D code. Since this example 
is much more demanding, we have found it important in 
order to increase the accuracy of our runs to perform a 
single initial maximal solve to reduce the initial gauge 
dynamics. The symmetries of this example are now not 
consistent with the evolution of just one octant. However, 
we still have reflection symmetry on the z = plane, so 
we evolve only the positive z half of the domain. The grid 
used in this case has 199 2 x 100 points and Ax = 0.4. 
The gauge conditions work well even in the presence of 
rotation: the shift drives the evolution to an almost static 
state as the system itself settles down to a Kerr BH. The 
metric functions (not shown) evolve in a similar way to 
those shown before, essentially freezing at late times. In 
Fig. 0, we show the extracted waveforms, now computed 
using the imaginary part of the Newman-Penrose quan- 
tity ^4 (e.g. H), which includes contributions from all 
i— modes at the same time. The results from the 2D and 
3D codes agree very closely, except for a slight phase shift 
due to slicing differences, until the 2D code becomes in- 
accurate and crashes. The 3D simulation continues well 
beyond this point, and is terminated at t = 120Af. 

Conclusions. We have extended recently developed 3D 
BH excision techniques, using a new class of live gauge 
conditions that dynamically drive the coordinates to a 
frame where the metric looks essentially static at late 
times, when the system itself settles to a stationary Kerr 
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FIG. 6. The solid line shows the £ = 2, m — waveform 
extracted at a radius of 5.45M for the even-parity distorted 
BH described in the text, while the dashed line shows the 
result of the same simulation carried out in the 2D code. We 
also show a fit to the two lowest QNM's of the BH for 2D and 
3D separately, using numerical data from t = 9M to t = 80M. 



BH. Our techniques have been tested on highly distorted, 
rotating BHs, and are shown to be very robust. For the 
first time, excision is tested with wave extraction, and 
waveforms are presented and verified. The results are 
shown to be more accurate, and much longer lived, than 
previous 3D simulations and even better resolved 2D sim- 
ulations of the same data. Such improvements in BH ex- 
cision are badly needed for more astrophysically realistic 
BH collision simulations, which are in progress and will 
be reported elsewhere. 

Acknowledgements. This work was supported by AEI. 
Calculations were performed using the Cactus code at 
AEI, NCSA, PSC, and RZG. 



[I] L. Smarr and J. York, Phys. Rev. D 17, 2529 (1978). 
[2] M. Alcubierre et al, (2000), gr-qc/0012079. 
[3] J. Baker et al, (2001), gr-qc/0102037. 
[4] R. Gomez et al, Phys. Rev. Lett. 80, 3915 (1998). 
[5] J. Thornburg, Class. Quan. Grav. 4, 1119 (1987). 
[6] E. Seidel and W.-M. Suen, Phys. Rev. Lett. 69, 1845 
(1992). 

[7] G. B. Cook et al, Phys. Rev. Lett 80, 2512 (1998). 

[8] M. Alcubierre and B. Briigmann, Phys. Rev. D 63, 

104006 (2001). 
[9] S. Brandt et al, Phys. Rev. Lett. 85, 5496 (2000). 
[10] L. E. Kidder, M. A. Scheel, and S. A. Teukolsky, (2001), 



0.00005 - 




2D(300x59) 

3L>(W5 J x0.5) 



FIG. 7. The solid line shows the imaginary part of tf>4 com- 
puted at r = 3.94M and = <j> — ir/4 for a rotating distorted 
BH obtained from our 3D code with excision, while the dash 
line shows the same quantity computed using a 2D code (this 
simulation crashes at t ~ 60M) . 



gr-qc/0105031. 
[11] P. Anninos et al, Phys. Rev. D 52, 2059 (1995). 
[12] G. E. Daues, Ph.D. thesis, Washington University, St. 

Louis, Missouri, 1996. 
[13] S. Brandt and E. Seidel, Phys. Rev. D 54, 1403 (1996). 
[14] S. Brandt and E. Seidel, Phys. Rev. D 52, 870 (1995). 
[15] K. Camarda and E. Seidel, Phys. Rev. D 59, 064026 

(1999), gr-qc/9805099. 
[16] J. Baker et al, Phys. Rev. D 62, 127701 (2000), gr- 

qc/9911017. 

[17] S. Brandt, K. Camarda, and E. Seidel, in Proceedings of 
the 8th Marcel Grossmann Meeting on General Relativity, 
edited by T. Piran (World Scientific, Singapore, 1999), 
pp. 741-743. 

[18] S. Brandt, K. Camarda, and E. Seidel, in preparation 
(unpublished) . 

[19] R. A. Matzner, M. F. Huq, and D. Shoemaker, Phys. 

Rev. D 59, 024015 (1999). 
[20] M. Shibata and T. Nakamura, Phys. Rev. D 52, 5428 

(1995) . 

[21] T. W. Baumgarte and S. L. Shapiro, Physical Review D 

59, 024007 (1999). 
[22] M. Alcubierre et al, Phys. Rev. D 61, 041501 (2000). 
[23] M. Alcubierre et al, Phys. Rev. D 62, 124011 (2000). 
[24] M. Alcubierre et al, Phys. Rev. D 62, 044034 (2000). 
[25] C. Bona, J. Masso, E. Seidel, and J. Stela, Phys. Rev. 

Lett. 75, 600 (1995). 
[26] J. Balakrishna et al, Class. Quant. Grav. 13, L135 

(1996) . 

[27] B. Briigmann, Int. J. Mod. Phys. D 8, 85 (1999). 

[28] M. Alcubierre et al, Class. Quant. Grav. 17, 2159 (2000). 



•5 



